//! Jacobi algorithm
use crate::mods::calculator::modulo_block;

/// For a small (usually prime, but not necessarily) prime p and >= 3 number,
/// compute Jacobi(p,num), which is -1, 0 or +1, using the following rules:
/// - Jacobi(x, y) = Jacobi(x mod y, y)
/// - Jacobi(0, y) = 0
/// - Jacobi(1, y) = 1
/// - Jacobi(2, y) = 0 if y is even, +1 if y is +/-1 mod 8, -1 if y = +/-3 mod 8
/// - Jacobi(x1*x2, y) = Jacobi(x1, y) * Jacobi(x2, y) (used with x1 = 2 & x1 = 4)
///
/// If x and y are both odd, then
/// > Jacobi(x, y) = Jacobi(y, x) * (-1 if x = y = 3 mod 4, +1 otherwise)
#[inline]
pub fn jacobi(mut p: i64, num: &[u64]) -> i64 {
    if p == 1 {
        1
    } else if p == 0 {
        0
    } else {
        let mut j: i64 = 1;
        let mut u: u64 = num[0];
        if (u & 1) == 0 {
            return 0;
        }
        /* Make p positive */
        if p < 0 {
            p = -p;
            let n8 = u & 7;
            if (n8 == 3) || (n8 == 7) {
                j = -j; /* 3 (11) mod 4 */
            }
        }

        /* First, get rid of factors of 2 in p */
        while (p & 3) == 0 {
            p >>= 2;
        }
        if (p & 1) == 0 {
            p >>= 1;
            if ((u ^ (u >> 1)) & 2) != 0 {
                j = -j; /* 3 (011) or 5 (101) mod 8 */
            }
        }
        if p == 1 {
            return j;
        }

        let mut p = p as u64;
        /* Then, apply quadratic reciprocity */
        if (p & u & 2) != 0 {
            j = -j;
        } /* p = u = 3 (mod 4)? */
        /* And reduce u mod p */
        u = modulo_block(num, p);

        /* Now compute Jacobi(u,p), u < p */
        while u != 0 {
            while (u & 3) == 0 {
                u >>= 2;
            }
            if (u & 1) == 0 {
                u >>= 1;
                if ((p ^ (p >> 1)) & 2) != 0 {
                    j = -j; /* 3 (011) or 5 (101) mod 8 */
                }
            }
            if u == 1 {
                return j;
            }
            /* Now both u and p are odd, so use quadratic reciprocity */
            if u < (p as u64) {
                let tmp = u;
                u = p;
                p = tmp;
                if (u & p & 2) != 0 {
                    j = -j;
                } /* u = p = 3 (mod 4)? */
            }
            /* Now u >= p, so it can be reduced */
            u %= p;
        }
        0
    }
}
